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Simulating frustrated quantum magnets is among the most challenging tasks in computational 
physics. We apply String-Bond States, a recently introduced ansatz which combines Tensor Net- 
works with Monte Carlo based methods, to the simulation of frustrated quantum systems in both 
two and three dimensions. We compare our results with existing results for unfrustrated and two- 
dimensional systems with open boundary conditions, and demonstrate that the method applies 
equally well to the simulation of frustrated systems with periodic boundaries in both two and three 
dimensions. 
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I. INTRODUCTION 

The simulation of correlated quantum spin systems is 
one of the central problems in condensed matter physics. 
The lack of exact solutions and the exponentially growing 
Hilbert space dimension motivate the need for numeri- 
cal methods for the simulation of such systems. Dur- 
ing the last decades, Quantum Monte Carlo (QMC) [l| 
and the Density Matrix Renormalization Group (DMRG) 
method 0, Q have arguably been the most successful 
methods for the accurate simulation of large quantum 
spin systems. Despite their huge success, both methods 
also have their limitations: The DMRG method gives 
extremely accurate results for one-dimensional (ID) sys- 
tems, but fails to simulate 2D systems similarly well; on 
the other hand, QMC can deal efficiently with 2D and 3D 
systems, but fails on frustrated (fermionic) quantum sys- 
tems due to the so-called "sign problem" . As frustrated 
quantum systems in two and three dimensions underlie 
some of the most interesting phenomena in condensed 
matter physics, methods which promise to overcome the 
previously mentioned limitations are of high interest. 

The natural generalization of the Matrix Product State 
(MPS) ansatz underlying DMRG to higher dimensional 
systems is given by Projected Entangled Pair States 
(PEPS) PEPS-based algorithms have been applied 

successfully, e.g., to the simulation of frustrated quantum 
spin systems or hardcore bosons in two dimensions [6|- 
LLOj . Yet, due to the scaling of resources the method is 
bound to two-dimensional systems with open boundaries, 
motivating the search for different tensor-network based 
algorithms pj]-[l8j . Recently, it has been proposed to use 
Monte Carlo sampling to enhance the possibilities of ten- 
sor network based methods, both in ID for DMRG [HI 
and, for appropriately chosen ansatz classes, in two and 
higher dimensions [20|, and their applicability to two- 
dimensional systems has been demonstrated |2CH22|. 

The String-Bond States (SBS) ansatz proposed in 
Ref. 20] generalizes the MPS ansatz to two and higher 



dimensions in a way which allows to employ Monte 
Carlo sampling to efficiently compute expectation values. 
While the MPS ansatz is inherently one-dimensional, 
SBS generalize it to higher dimensional lattices by plac- 
ing several one-dimensional structures atop of each other, 
e.g., along the axes, the diagonals, and in loops between 
adjacent neighbors, thus allowing for arbitrary correla- 
tions between any group of spins without sacrificing the 
advantages of the one-dimensional structure. 

In this paper, we demonstrate the applicability of 
the SBS ansatz to the simulations of two- and three- 
dimensional frustrated quantum systems. In two dimen- 
sions, we apply it to the simulation of the frustrated J±- J 2 
model, where we find that for open boundary condition 
(OBC), SBS reproduce well both the energies and the 
structure of correlations obtained using the general PEPS 
ansatz. Moreover, SBS also allow us to simulate systems 
with periodic boundaries (PBC) with similar accuracy, 
and we find that the behavior of the low-energy regime 
of the system in the transition region J2/J1 ~ 0.6 (chang- 
ing from Neel to columnar order) differs significantly for 
OBC and PBC. 

Second, we apply the SBS ansatz to the simulation of 
3D frustrated spin systems. To benchmark the ansatz, we 
compare results for the 3D Ising model with transverse 
field to results obtained using QMC. Then, we apply it 
to the simulation of a three-dimensional frustrated quan- 
tum spin system on up to 6 x 6 x 6 = 216 qubits, where 
we observe a performance comparable to that in two di- 
mensions. This demonstrates the ability of the method 
to simulate frustrated quantum spin systems in both two 
and three dimension and with periodic boundaries. 



II. THE STRING-BOND STATE ANSATZ 

String-bond states have been proposed as a variatonal 
class of states for which expectation values of local ob- 
servables can be computed efficiently using Monte Carlo 
sampling. Monte Carlo sampling allows to compute an 



expectation value ^2p(n)f(n) over a probability distri- 
bution p by generating a sample {ni, 77,2, . . .} drawn from 
p(n) and averaging / over this sample. Now for any ob- 
servable O, we can rewrite its expectation value in a state 

1*0) as 



<V>|0|V> = X>l»>MOM = £>(n) 



(n|V) 



(1) 



where p(n) = Kn^)! 2 and \n) is an orthonormal basis. 
Thus, whenever (n\%j)) and (n\0\ip) can be computed effi- 
ciently, (-01 01-0) can be evaluated efficiently using Monte 
Carlo sampling. 

We are interested in investigating systems consisting of 
TV spins with a Hilbert space (C d )® N , and we thus choose 
the basis \n) = \n±, . . . , un) to be a product (i.e. local) 
basis of the system. In order to do efficient Monte Carlo 
sampling in this basis, we need that (n\i/j) and (n\0\tp) 
can be computed efficiently. The second requirement can 
be reduced to computing a few overlaps (h\ip) whenever 
O = ^DjzPjz with Dk diagonal, permutations, and 
the set of fc's sufficiently small, since then (n\0\ip) = 
^2 k Dk(n)(nk\^)^ with (n\k = {n\Pk another local basis 
state. In particular, this holds for local O (where local 
means small support, as e.g. the terms in a Hamiltonian 
or two-point correlation functions) and tensor products 
of Paulis, for instance the Jordan- Wigner transform of 
fermionic hopping terms, or string order parameters. 

Thus, in order to be able to apply Monte Carlo sam- 
pling, we need to find classes of states for which the 
overlap {n\ijj) can be computed efficiently. We choose 
(n\i/j) to be a product of efficiently computable func- 
tions f s (with s = l,...,*?) defined on subsets M s C 
{!,..., N} of spins, 



(2) 



Here, nj^ s contains the state of all spins in the subset 
M s . Note that the subsets M s should be overlapping as 
otherwise they just describe a product state. 

Our choice of the f s will be such as to generalize Ma- 
trix Product States (MPS) to higher dimensional sys- 
tems. An MPS of bond dimension D is given by 

|V> = Yl tr [ M ni ' ' ' M n N ] \ni,..., n N ) (3) 



where My are D x D matrices. In order to generalize 
MPS in the spirit of the ansatz (j2j), we choose each f s 
such that 



f s (n il ,...,n il ) = tr 



s.l 



s,l 



(4) 



to be a trace of matrix products. Here, zi, . . . , i\ denotes 
the spins in the corresponding subset M 3 \ note that this 
imposes an ordering on these sets. Clearly, this definition 
includes MPS themselves, since we can choose only one 
Af s = {l,...,N}. 

In defining SBS on higher dimensional systems, the 
choice of the subsets M s (called "strings" furtheron, as 
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FIG. 1: (Color online). String patterns used in the simula- 
tions, a) The basic lines pattern. It can be enhanced by the 
b) diagonals pattern and by the c) loops pattern, which help 
to improve the control over diagonal and four-body correla- 
tions, respectively. 



they impose a one-dimensional ordering in the spirit of 
MPS) is of central importance. The idea is that the string 
pattern should reflect the geometry of the system in such 
a way that spins which are closely coupled by the Hamil- 
tonian are rather closely connected by a string. For a 2D 
square lattice, a natural choice is to first put strings on 
all rows (i.e., one row forms one string, corresponding to 
a product of MPS on rows) and then connect the rows by 
additionally placing one string per column. We call this 
pattern, as illustrated in Fig. QJi, lines. The lines pattern 
can be enhanced in two different ways by putting addi- 
tional strings: First, one can put strings on all diagonals 
of the lattice (Fig.HJ)), and second, one can choose strings 
which form small loops, encompassing all elementary pla- 
quettes (i.e., blocks of 2 x 2 spins, Fig. [It; cf. [21] for a 
generalization of this ansatz); both of these extensions al- 
low for a better control of the correlations with diagonal 
neighbors. The patterns generalize straightforwardly to 
lattices in 3D or with different geometries. Note that by 
continuously adding strings, we will eventually be able 
to describe all states as SBS, as can be seen by putting 
one long snail-like string on the lattice (i.e., describing 
the whole state as an MPS). Clearly, for good practical 
results, the strings should be chosen such that the rel- 
evant states are well approximated at an early stage of 
the pattern. 

The computational resources of SBS scale favorable as 
compared to PEPS: For each string, a matrix trace (|4]) 
has to be computed which takes resources ID 3 (ID 2 for 
OBC), with / < TV the length of the string. This has to 
be multiplied by the number of strings *S, giving a com- 
putational cost of 0(SND 3 ). In particular, the scaling 
in the accuracy parameter D compares favorably to the 
D 10 (D 18 ) scaling of the PEPS method for OBC (PBC). 

Let us briefly note that although we motivated SBS 
as a higher-dimensional generalization of MPS, one can 
also regard them as a specialized case of PEPS. PEPS 
form the most natural generalization of MPS to two di- 
mensions they are known to approximate the states 
of interest well [23|, [24| , and have been applied success- 
fully in numerical simulations [H, 0]- However, the scal- 
ing in the accuracy parameter is rather bad, preventing 
the application of PEPS to problems beyond 2D systems 
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with OBC (note, however, that iPEPS have been applied 
successfully to investigate 2D systems in the thermody- 
namic limit 0, [Iq|). One way to resolve this problem 
is to look for subclasses of PEPS which allow for more 
efficient algorithms. Indeed, SBS form such a subclass of 
PEPS [U: While general PEPS are described by tensor 
networks with general tensors T^^, SBS with a lines 
pattern have tensors of the form Ai a pBi 7 §. Note, how- 
ever, that the structure of the tensors gets more and more 
rich as one places additional strings on the lattice, and 
thus, SBS can only be embedded in PEPS at a cost ex- 
ponential in the number of strings; moreover, since SBS 
computations scale much more favorably in the accuracy 
parameter, even for a basic lines pattern SBS can out- 
perform PEPS as they can reach much larger D's. 



III. VARIATIONAL METHOD USING 
STRING-BOND STATES 

In the previous section, we have introduced string- 
bond states (SBS) as a class of states which generalize 
MPS to two- and higher dimensional systems while allow- 
ing for an efficient computation of expectation values. In 
this section, we will show how SBS can be used to build a 
variational algorithm for simulating the ground states of 
quantum spin systems. Although the ability to efficiently 
compute expectation values is a necessary criterion, it is 
not sufficient: One also needs an efficient and practical 
way to evolve the SBS towards the ground state. 

The basic idea of a variational algorithm based on SBS 
is to fix a family of SBS (i.e., fix a certain string pattern 
and the dimension D of the underlying matrices) and try 
to find the state within this family which minimizes the 
energy of a given local Hamiltonian. Similar to DMRG 
or the variational method over PEPS, we will carry out 
the optimization in a local fashion: We start from some 
SBS, described by a number of three- index tensors M as 
in (HJ), select one of the tensors - let us call it A - and try 
to minimize the energy with respect to this tensor while 
keeping the others fixed. This procedure is repeated for 
all tensors over and over until the energy converges, i.e. 
a minimum within the family of states is reached. 

To determine how to change the selected tensor A such 
as to minimize the energy, we use the linearity of the 
string-bond states in the tensor A to be optimized, 



(ipA\H\iP A ) (A\X\A) 



(iPa\iPa) 



(A\Y\A) 



(5) 



where we have explicitly denoted the dependence of 
the string-bond state \?p A } on A. (A\X\A) denotes a 
quadratic form in A, where \A) is the vectorized form of 
A, i.e. A(ij k ) = A%j, and we use boldface to avoid con- 
fusion with vectors in state space. Minimizing ((5]) with 
respect to A is a generalized eigenvalue problem and can 
be solved efficiently. 

In order to sample X and Y, define vectors \a n ) and 



\b n ) via the linear functionals 



(n\tpA) 

( n \i>A ) ' 



(6) 



where Aq is the initial value of the tensor A. It follows 
that the matrices X and Y in §5§ can be expressed as 



X = 



(n)\b n ) (a n \ , Y = ^Po(n)\b n )(b n \ , (7) 



where po(n) oc K^I^Ao)! 2 ? an d thus determined by Monte 
Carlo sampling of and |6 ri )(6 Tl |, respectively. 

Note that by virtue of this definition, we obtain the nor- 
malization (AqII^IAq) = 1. 

However, there is a major problem with the approach 
of solving the generalized eigenvalue problem: Monte 
Carlo sampling X and Y is relatively inaccurate as com- 
pared to e.g. the approximate contraction as done in the 
PEPS algorithm [4], and moreover, our estimates of X 
and Y get less and less accurate for A's far away from 
Ao as we have sampled with respect to the distribu- 
tion at A = Aq. Specifically, already small errors in Y 
might lead to completely wrong minima for the general- 
ized eigenvalue problem. While for the PEPS algorithm, 
this problem can be successfully overcome by truncating 
small eigenvalues of Y , this is impractical for a method 
based on Monte Carlo due to the comparatively large 
error. 

To overcome this problem, we do not solve the gen- 
eralized eigenvalue problem to compute the new A, but 
rather compute the gradient of the energy with respect 
to A and change A slightly along this gradient such as to 
decrease the energy. First, this accounts for the fact that 
our sample of X and Y , Eq. (0, is most accurate around 
Ao, and as we will see, it moreover yields a formula where 
neither X nor Y appear in the denominator, such that 
small absolute errors remain small. Another advantage 
will be that it is possible to gain a considerable speed- 
up when sampling all the gradients simultaneously and 
change all the tensors along their gradient simultaneously 
- this is possible since the gradients decouple to first or- 
der. 

In order to determine the gradient of the energy with 
respect to A around Ao, consider a small variation A = 
A + eB (with e«l): 



Ao+eB 



(A +6B\X\Ao + 6B) 

{Ao + eB\Y\A + eB) 

(A \X\Aq) + 2e Re j(B\X\A )} + 0(e 2 ) 

l + 2eRe[(B|r|A )] + O(e 2 ) 
E^ Ao ) + 2eRe[(B\X\A )} 

- 2e(A \X\A )Re [(B\Y\A )] + 0(e 2 ) 



where we have used the normalization (Ao|l^|Ao) 
Thus, the gradient turns out to be 



1. 



V A E(^ A )\ =2[X\A )-E^ Ao )Y\A )] (8) 



[using (A |X|A ) = E^Aq)]- Substituting the sampling 
formulas ([7]) for X and Y and using that (b n \A ) = 1 
[Eq. flg)], we finally obtain that 

V a E^ a )\ a ^ Aq = 2^> (n)|6 ri ) [E n - E(^ Ao )] , 



where we have defined 



E n := (a n \A ) 



and the energy can be computed as 



While the previous derivation holds for any ansatz 
where \iP)a is linear in A, there are some additional tricks 
which can be applied in the case of SBS to save compu- 
tation time. To this end, note that all we have to know 
are |6 n ), E ni and the ratio po(n)/po(m) (this is sufficient 
to generate a random walk) . For a particular tensor ■ , 
the dependence of {ti^a) on A (where ua denotes the 
state of the spin associated with A) can be expressed as 

(n\1> A ) = ti[A nA X(n)]c(n) 

where X(n) is the product of all other matrices on the 
string containing A as a function of the state n of the 
spins, and c(n) contains the contributions from all other 
strings. Thus, we have that 



(b n \A) 



{n\^ A ) _ ti[A nA X(n)]c(n) 



(n\^Ao) tr[A^X(n)]c(n) 
and therefore (with s the physical spin index) 



t X(n)t 



trK-X(n)] ' 



This means that in order to compute |6 n ) for a given ten- 
sor A, one only has to consider the string which contains 
A, instead of having to look at all the strings. 

Similarly, in order to compute E ni one can exploit 
that for local Hamiltonians, string-order operators, etc., 
(n\H = ^2 meM f(m)(m\ where A4 has only few ele- 
ments, and e.g. for local Hamiltonians on a 2D lattice, 
each m G M only differs at two adjacent sites from n. 
Thus, 



E n = 



(n\H\iP Ao ) 
(n\ipA ) 



y, (m\ip Ao ) 



can again be computed as the ratio of the matrix prod- 
uct traces for only the two strings on which m and n 
differ, again reducing the computational effort. Comput- 
ing po(m) / po(n) also allows for optimizations, depending 
on the way the new configuration m is constructed start- 
ing from n. For the simplest scenario where only a single 



spin is flipped, again only the strings containing this very 
spin have to be considered, and similarly if e.g. a pair of 
spins is being flipped. 

Finally, we gain a speed-up by computing the gradients 
for all tensors simultaneously; this is reasonable since 
the joint gradient of all tensors is nothing but the direct 
sum of the individual gradients, thus, changing all the 
tensor in direction opposite to the gradient by a small 
amount will decrease the energy to leading order. In this 
case, computation time is saved by the fact that the same 
sample drawn from po(n) can be used, and that E n has 
to be computed only once. 

The full algorithms looks as follows: Fix a string pat- 
tern and corresponding bond dimensions, and choose ini- 
tial configurations for all tensors. Then, iterate the fol- 
lowing: 1) Compute the energy and its gradient with 
respect to all tensors. 2) Change all the tensors by some 
small amount in the direction given by the gradient. 3) 
Start over at 1) with the modified tensors. Iterate this 
until the change in energy becomes smaller than some 
threshold and declare convergence. In order to ensure 
that the step along the gradient is small enough, it is 
advisable to normalize the gradients such that the step 
remains small even for steep gradients. 

Instead of declaring convergence of the algorithm when 
the energy does not change any more, one can try to 
increase the precision and see whether this leads to a 
further improvement in energy, and only declare conver- 
gence if it doesn't. There are three possiblities to do so: 
First, one can increase the length of the Monte Carlo 
sample used for computing the energy and the gradient, 
second, one can try to decrease the stepwidth used to up- 
date the tensors along the gradient, and finally, one can 
try to extend the variational family of states either by in- 
creasing the bond dimension or by adding extra strings. 
In all cases, it is advisable to use the previously obtained 
optimum as the initial state. 



IV. NUMERICAL RESULTS 

In the following, we present numerical results obtained 
for two- and three-dimensional frustrated spin systems 
using string-bond states. In all cases, the Monte Carlo 
sampling was carried out using single spin flip, or adja- 
cent spin swap, Metropolis updates. The autocorrelation 
time was at most 100 updates (for the structure factor of 
the J1-J2 model in the frustrated regime), and consider- 
ably less for local observables or non- frustrated models, 
even in 3D. This allowed us to choose the Monte Carlo 
samples sufficiently long such that in all cases, the error 
bars were below what could be illustrated in the plots 
(local observables to at least 0.1%, and non-local observ- 
ables to at least 1% accuracy). Note, however, that this 
only means that we have good control over the error we 
make in measuring observables on the given variational 
state; the major (and not so well controlled) error source 
in the method is thus the ability of the ansatz class to 
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FIG. 2: (Color online.) Relative error in ground state energy 
as a function of the string pattern for a the J\-Ji model on a 
6x6 OBC lattice. From left to right: lines with D = 2, 4, 6, 8 
(labelled L2, . . . , Ls), Zmes with D = 8 together with diago- 
nals with D = 2, 4, 6, 8 (labelled D2, • • • , Ds), and finally lines 
with D = 8, diagonals with D = 8, and loops (D = 4). The 
lower (red) points are for the Heisenberg model, J2/J1 = 0, 
and the upper (red) blue for J 2 /Ji = l. In the highly frus- 
trated regime, adding diagonals and loops leads to a signifi- 
cant improvement. Note that an extrapolation is difficult to 
perform, as it is unclear how the accuracy will scale in the 
string pattern. 



correctly describe the ground state, together with the 
question as to whether the variational method converges 
to the optimal state within the class. 



A. Simulation in 2D: The J\ — J2 model 

We have applied SBS to the simulation of the so-called 
J 1- J 2 model, 



H 



JU2 



<i,j> 



J2 sr^ 



where < z, j > denotes nearest neighbors in a 2D square 
lattice, and <C i, j ^> nearest neighbors along the diago- 
nal. This model arises e.g. in the context of the Hubbard 
model which is believed to underly high-temperature su- 
perconductivity [IK , and has become one of the paradig- 
matic models to understand quantum phase transitions 
in frustrated spin systems [26j |. 

For the simulation, we started from the patterns lines, 
then added diagonals, and finally loops. Fig. [2] shows 
how the energy improves as D is increased and additional 
strings are added, for J2/J1 = 1; as on one can see, the 
improvement due to additional strings depends on the 
model under consideration. Note that for our simula- 
tions, we have used the SU(2) invariance of the model, 
which implies that we can project our ansatz into the 
spin subspace (as there is a ground state with spin 0). 
This can be understood as an SBS with one additional 
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FIG. 3: (Color online). Energy comparison for the J\-Ji 
model for OBC lattices of size 4x4 (compared to the exact 
energies), 6x6, and 10 x 10 (both compared to the PEPS 
energies Q). 
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FIG. 4: (Color online.) Energy comparison for the J1-J2 
model for PBC lattices of size 4 x 4, 6 x 6 (both compared to 
exact energies [26]), and 10 x 10 (where there is no data to 
compare with). 



string which covers the whole lattice and enforces S z = 0. 
In practice, we achieve the restriction by sampling from 
the S z = subspace: we start from a configuration in 
this subspace and create new configurations by swapping 
a randomly chosen pair of spins; we have observed that 
this restriction led to a significant improvement in energy. 

In Fig. [3j we show results for the ground-state energy 
of the J\- J2 model on lattices of size 4 x 4, 6 x 6, and 
10 x 10 with open boundaries, which we compare with 
the values obtained using exact diagonalization (4 x 4) 
and the PEPS method [7] (6 x 6, 10 x 10). Fig.Hshows 
the same numbers for the case of periodic boundaries, 
compared to the exact numbers [26| (4x4, 6x6). Let 



J 2 IJ X = 

\S(4> 

PEPS OBC 



SBS OBC 
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FIG. 5: (Color online). The absolute value of the structure factor S(4> x , <p y ) as defined in Eq. Q computed for the J1-J2 model 
on a 10 x 10 lattice as a function of the ratio J2/J1. The plot compares the results obtained on an OBC lattice using PEPS 
with both the OBC and the PBC result found using SBS. One finds that for OBC, SBS reproduce the characteristics of the 
PEPS results, and there is the signature of a intermediate glassy phase around J2/J1 — 0.6. For PBC, on the contrary, there is 
no signature of an intermediate phase, which is missing for PBC. Note that this observation should be taken with care, as the 
SBS energies are typically a few percent above the PEPS. The wave- like artifacts which can be seen especially for OBC around 
J2/J1 = 0.6 are probably due to the fact that the string pattern has preferred axes. 



us note that for 10 x 10 lattices, there are no numbers 
available to compare with. Typical D's were Aine = 
Aliag = 6 for up to 6 x 6 and 8 to 10 for 10 x 10, and 



4. 



The relative errors in energy corresponding to Figs. [3] 
and [H are shown in Fig. [6l While the energies obtained 
using SBS are above the exact/PEPS data, the error does 
not seem to depend on the system size or the choice of 
boundaries, which suggests that the method should be 
equally applicable to larger and PBC systems. 

Let us now see whether SBS can reproduce the corre- 
lation functions of the J1-J2 model. To this end, we use 
the structure factor 



(9) 



S(<p) is the Fourier transform of the two-point correla- 



tion functions (ovr^m)? i- e -? it reveals information about 
the relative alig nment of the spins, this is, the order of 
the system [27f. The results for for PEPS with OBC, 
SBS with OBC, and SBS with PBC is diplayed in Fig. El 
Note that the OBC results exhibit the same character- 
istic properties for both PEPS and SBS, while the SBS 
results for PBC are significantly different in the region 
around J2/J1 = 0.6, where the behavior of the model is 
not yet fully understood. It is believed that in this re- 
gion, the system is in some kind of glassy phase. While 
a signature of this phase can be seen in the case of OBC 
both for the PEPS and the SBS data, the same signature 
is completely absent in the case of periodic boundaries. 
While this seems to suggest that the behavior of the sys- 
tem in that region might be different for OBC and PBC, 
we would like to stress that this is obtained from config- 
urations with energies clearly above the exact and PEPS 
data, and thus should be treated with care. 
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FIG. 6: (Color online). Relative error for the comparisons in 
Figures [3] and [4] As one can see, the accuracy stays constant 
when increasing the lattice size. 



B. Three-dimensional systems 

While some variational methods based on tensor net- 
works such as PEPS 0, MERA [l8j, or Monte-Carlo 
based ansatzes such as SBS [2(j[ or EPS (2l| have been 
shown to be able to simulate two-dimensional frustrated 
quantum systems, none of the previous methods has yet 
be applied to the simulation of systems in three dimen- 
sions. In the following, we give results of SBS simulations 
for three-dimensional frustrated quantum systems which 
are comparable to those obtained in two dimensions. 

The 3D simulations are based on the lines pattern on 
a 3D lattice with PBC, and D = 6. To benchmark 
the method, we have simulated the 3D Ising model with 
transverse field, H = ^ ZiZj + B ^ Xi, on an 8 x 8 x 8 
PBC lattice, and compared the result to QMC Simula- 
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FIG. 7: (Color online). Results for the 3D Ising model in 
transverse field on an 8 x 8 x 8 PBC lattice. The plot shows 
the magnetization obtained with SBS using the lines pattern, 
compared to QMC data obtained using ALPS [28, [29| , and to 
the mean field solution. The inset shows the relative error in 
the energy. 



FIG. 8: (Color online). Squared magnetization (M%) for the 
3D transverse Ising model on an 8 x 8 x 8 PBC lattice, com- 
paring data obtained with SBS using the lines pattern, QMC 
data obtained using ALPS [28l [29j], and the mean field so- 
lution. The inset shows a log-log-plot close to the critical 
point. See text for a comment of the fluctuations which can 
be observed around the critical point. 



tions carried out using the ALPS package [28|, [29j, as 
well as mean field data. Fig. [7] shows the magnetization 
along x and the relative error in energy (inset) as a func- 
tion of the field B, and Fig. [8] the magnetization squared 
along the Ising coupling, (M%) = Y,ij( z i z j) / N<2 • Note 
that the method becomes unstable close to the critical 
point and frequently gives too large values for (M|). 
This is not a problem of the Monte Carlo sampling, 
which yields (M^) with an accuracy of about 1% (with 
1.6 x 10 6 Metropolis updates, where (M%) is sampled on 
every 100th configuration). Rather, this effect is due to 
the fact that variational methods using MPS and related 
ansatzed such as SBS generally tend to break symme- 
try close to the critical point even in ID, as has been 
also observed elsewhere [30j. This can be understood 
in two ways: Firstly, the entanglement entropy of the 
ground state diverges at the critical point, so that the 
ground state cannot be exactly reproduced by states such 
as MPS or SBS which obey an area law, thus driving 
the ansatz into symmetry-broken solutions with slightly 
higher energy but less entanglement. Secondly, varia- 
tional ansatzes have a general tendency to break symme- 
tries as this corresponds to having less (connected) long- 
range correlations, and establishing such correlations is 
difficult to accomplish by doing local optimizations. E.g., 
in the most extreme case, once the matrices in the MPS 
or SBS do not have full rank any more, the subspace 
not used by the matrix is lost for the optimization as 
it cannot be seen any more by local variations, and in 
particular by a gradient search. 

After having tested our 3D algorithm on the Ising 
model, we have subsequently applied SBS to simulate a 
frustrated XX model in a transverse field on a 3D square 
lattice, 



where <i,j> denotes nearest neighbors on the 3D square, 
and Jij = ±1 is chosen such that the system is frustrated 
around every plaquette, as illustrated in Fig. [9] There are 
several reasons for chooosing this model: First, it is frus- 
trated and thus cannot be simulated by QMC due to the 
sign problem. Second, its lower symmetry as compared 
to an SU(2) invariant model makes it easier to simulate. 
Finally, for this model, the z magnetization S z is a good 
quantum number. Thus, the behavior of the model can 
be completely understood if the minimal energy E m for 
fixed m = S z at zero field is known: The minimal en- 
ergy within each subspace with given magnetization m 
decreases linearly with the field, E rn (B) = Em — Bm, 
and the magnetization m at a given B is the one for 




<i,j> i 



(10) 



FIG. 9: (Color online). Coupling pattern for the 3D frus- 
trated XX model Eq. ([TO]) , illustrated for an elementary cell 
of size 2x1x2. The thick red edges represent ferromag- 
netic couplings Jij = —1, while the other edges correspond 
to antiferromagnetic couplings, Jij = 1. Note that the role 
of ferromagnetic and antiferromagnetic couplings, as well as 
the axes of the model, can be swapped by local a z transfor- 
mations. 
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magnetic field 



FIG. 10: (Color online). Magnetization of the 3D frustrated 
model (|1Q|) as a function of the external field for a 6 x 6 x 6 
lattice, for SBS and mean field. The inset shows the improve- 
ment in energy of SBS relative to mean field. 



which E m (B) becomes minimal. 

We have computed the E m for the model ([TO]) on a 
6x6x6 lattice and from this data determined the ground 
state energy and the magnetization as a function of the 
field. The results are shown in Fig.[10j where we compare 
it to mean field data, which we also used to bootstrap the 
SBS ansatz. We found that most of the improvement is 
already obtained for D = 2 (D = 1 being mean field), 
and for D = 6, the method was fully converged. 

In order to estimate the performance of the ansatz, we 
have compared both mean field and SBS to the exact 
solution on a 2 x 2 x 4 lattice. The results are shown in 
Fig. [TTJ While both energy and magnetization are still 
away from the exact solution, the values obtained using 
SBS are significantly more accurate than the mean field 
solution. 




FIG. 11: (Color online). Benchmark for the 3D frustrated 
model (|10p on a 2 x 2 x 4 lattice: We compare the exact val- 
ues with data obtained using mean field and SBS. The figure 
shows the magnetization as a function of the magnetic field, 
and the inset the improvement in energy relative to mean 
field. 



V. CONCLUSIONS 

In this work, we have presented numerical results ob- 
tained with the recently introduced String-Bond State 
(SBS) ansatz for frustrated quantum spin systems in 
both two and three dimensions and for open and periodic 
boundaries. While the results obtained for 2D OBC sys- 
tems were above the results found using PEPS, the more 
favorable scaling of the method allowed us to go beyond 
2D and OBC and obtain similarly accurate results for 2D 
PBC and 3D frustrated systems, which often cannot be 
simulated otherwise. 

The computational resources needed for the simulation 
are moderate, as the contraction of the strings scales only 
with D 3 , and the Z}'s used are much smaller than those 
in DMRG; typical simulations for the J1-J2 model took 
less than two days using a MATLAB code on a single 
processor. The method allows for parallelization in eval- 
uating energy and gradient on the Monte Carlo sample 
with low interprocess communication. Optimizations are 
possible with respect to caching contracted strings and 
reusing them in consecutive Monte Carlo samples, as well 
as in reusing Monte Carlo samples after small updates. 

There are two main challenges in the implementation 
of the algorithm. First, one needs a systematic way 
of growing the string pattern which is suitable for the 
problem at hand. As one can see in Fig. [2j the same 
string patterns lead to different improvements depend- 
ing on the underlying model. Related to this, the per- 
formance of the method on non-SU(2) invariant models 
will also depend on the choice of the local basis in which 
the sampling is performed, since this will affect the prob- 
ability distribution sampled over. The second important 
point is to choose the proper initial state for the opti- 
mization. In particular, we have observed for the three- 
dimensional frustrated XX model presented in the paper 
that the algorithm performs much better when starting 
from the mean field solution as compared to a random 
initial state. Here, it seems that the important informa- 
tion is the proper sign of the wavefunction rather than 
the amplitude, as the latter can be easily changed by the 
gradient flow. (Note however that the performance for 
the J 1- J 2 model did not depend on the choice of the ini- 
tial state.) The proper choice of the sign pattern will 
likely also pose a central challenge when applying SBS 
to fermionic systems; note however that this might be 
overcome by using fermionic SBS, analogous to fermionic 
PEPS HH, |32|, instead of mapping the system to spins 
via a Jordan- Wigner transform. 
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